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Technology,  under  the  Defense  Research  Sciences  Program,  Project  2304, 
Mathematical  and  Information  Sciences,  Task  N1 ,  Computational  Aspects 
of  Fluid  and  Structural  Mechanics.  It  constitutes  the  Final  Report  of 
Work  Unit  2304N112,  Duct  Acoustics  Research,  describing  work  accomplished 
between  December  1975  and  May  1981. 
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SECTION  I 
INTRODUCTION 

Duct  acoustics  continues  to  be  of  great.  Interest  to  the  aircraft 
industry  because  of  the  need  of  obtaining  quieter  aircraft  by  inserting 
sound  absorbing  linings  In  jet  engines.  When  considering  a  reasonably 
detailed  representation  of  an  aircraft  engine,  one  obtains  a  system  of 
linear  differential  equations  with  nonconstant  coefficients  to  be  solved 
in  a  region  with  a  complicated  nonuniform  geometry.  Earlier,  we  used  a 
finite  difference  method  with  a  conformal  map  to  determine  sound  pressure 
levels  within  nonuniform  ducts  in  the  absence  of  flow  (Referer.ces  1  and  2). 
However,  for  all  but  a  few  duct  geometries  the  conformal  map  must  be 
determined  numerically,  and  this  can  be  a  complicated  problem  in  Itself. 

An  integral  equation  method  was  used  in  Reference  3  to  solve  the  no-flow 
equations  in  both  uniform  and  nonuniform  ducts.  This  method  is  very 
efficient  in  terms  of  computational  time  and  does  not  require  a  mapping 
to  a  uniform  geometry.  However,  for  the  problem  of  flow  in  a  nonuniform 
duct,  the  acoustic  equations  have  variable  coefficients.  The  fundamental 
solution,  which  is  needed  in  the  integral  equation  method,  is  not  known. 
With  the  finite  element  method,  we  are  able  to  handle  complicated  non- 
uniform  geometries  and  also  take  advantage  of  the  high  order  convergence 
of  the  higher  order  finite  element  method?*.  For  these  reasons,  we  have 
decided  to  evaluate  the  finite  element  method  for  solving  problems  in  duct 
acoustics.  Because  of  the  great  versatility  of  the  finite  element  method, 
there  are  many  alternatives  to  consider.  For  instance,  a  least  squares 
approach  could  be  used,  and  the  resulting  positive  definite  system  of 
equations  (called  the  stiffness  matrix  by  the  structural  engineers)  could 
be  solved  by  an  iterative  method.  The  use  of  an  iterative  method  had  been 
considered  and  rejected  for  finite  differences  because  the  finite  difference 
system  of  equations  was  not  positive  definite  and,  hence,  not  compatible 
with  an  iterative  solution  technique.  However,  the  system  of  equations 
which  arises  from  a  least  squares  finite  element  approach  is  positive 
definite  and  can  be  solved  iteratively. 
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A  problem  with  the  least  squares  approach  Is  that  it  introduces  an 
additional  source  of  error  which  does  not  occur  with  a  Galerkin  finite 
element  approach.  So  perhaps  a  Galerkin  approach  should  be  used.  This 
can  be  viewed  as  a  type  of  weak  solution  of  the  differential  equations 
In  question  or  as  a  method  of  weighted  residuals.  The  system  of  equations 
which  results  when  a  Galerkin  approach  is  used  to  solve  the  Helmholtz 
equation  is  not  positive  definite  and,  therefore,  cannot  be  solved 
Iteratively.  If  a  direct  solver  is  used,  the  decision  has  to  be  made 
whether  to  make  it  in-core  or  out-of-core  on  the  computer  at  hand. 

Once  the  basic  approach  has  been  decided,  a  choice  of  basis  or  shape 
functions  must  be  made.  The  two  choices  considered  in  this  paper  are 
piecewise  linear  and  piecewise  quadratic  polynomials.  The  accuracy  of 
these  two  different  systems  of  basis  functions  are  compared  as  well  as 
their  respective  computational  times.  Related  to  the  choice  of  basis 
function  is  the  way  that  the  duct  geometry  is  divided  into  elements. 

One  possibility  considered  was  the  use  of  rectangles  and  triangles  with  the 
triangles  used  to  improve  the  approximation  of  the  boundary  of  a  non- 
uniform  duct.  A  second  approach  consists  of  using  only  quadrilaterals 
with  perhaps  curved  sides  to  approximate  the  boundary.  Results  were 
consistent  with  others  (Reference  4)  in  the  discovery  that,  in  general, 
triangles  have  lower  accuracy  than  rectangles  or  quadrilaterals  for 
many  classes  of  problems.  Since  t.he  polynomials  and  products  of  poly¬ 
nomials  which  comprise  the  basis  functions  need  to  be  integrated,  It  has 
to  be  decided  whether  to  integrate  exactly  or  to  use  numerical  quadrature. 
These  questions  and  alternatives  have  been  addressed  and  will  be  discussed 
in  subsequent  sections  of  this  paper. 

One  of  the  most  difficult  decisions  to  be  made  when  solving  duct 
acoustics  problems  in  nonuniform  ducts  is  choosing  the  type  of  flow. 

The  flow  considered  Is  that  which  corresponds  to  uniform  mean  flow  in  a 
uniform  duct.  This  flow  field  then  provides  the  coefficients  in  the 
acoustic  equations. 
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As  the  results  Indicate,  the  finite  element  method  offers  simplicity 
of  use  once  the  Initial  programming  and  bookkeeping  efforts  are 
completed  and  Is  accurate  when  compared  to  competing  methods  such  as 
finite  differences.  The  advantage  that  nonuniform  geometries  are  more 
easily  handled  by  finite  elements  versus  finite  differences  Is  not  as 
great  when  the  mapping  function  Is  used  to  compute  the  flow  field.  This 
Is  because  the  same  mapping  function  could  be  used  to  map  the  nonuniform 
duct  to  a  rectangular  duct  with  the  acoustic  equations  on  the  rectangle 
then  solved  by  finite  differences. 

In  Section  II  the  derivation  of  the  acoustic  equations  and  boundary 
conditions  will  be  described.  In  Section  III  the  finite  element  solution 
method  will  be  discussed  and  In  Section  IV  the  results  will  be  presented. 
For  other  references  using  the  finite  element  method  to  solve  duct 
acoustics  problems,  see  References  5-8. 
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SECTION  II 

DERIVATION  OF  THE  EQUATIONS 
1 .  DIFFERENTIAL  EQUATIONS 

If  viscous  and  heat  transfer  terms  are  neglected ,  the  equations  of 
motion  for  an  ideal  gas  are: 

Momentum  p'  D($' )/Dt  ■  -Vp'  (a) 

Continuity  D(p')/Dt  +  p‘(V*$' )  »  0  (b)  (1) 

Energy  and  state  p'  *  const  pa  (c) 

where  p'  Is  the  density,  1}'  Is  the  velocity  vector,  and  p'  Is  the  pressure. 


Equations  1  can  be  combined  to  form  the  usual  conservation  of  mass 
and  momentum  equations: 

If  +  *'-vp'  +  *  0 

p'  <|f  ♦  (2) 

Because  of  the  application  to  aircraft  engines,  these  equations  have  been 
written  in  cylindrical  coordinates  and  become: 


at 


“•If 


+  V 


3r 


vT 

r 


pu1 

-  ♦  av 

+  I  pL  ■ 

.  vll  . 

n 

L82 

3r 

r  36 

V 

+  u' 

|yi  +  v 

,9u'  .  vf 

'  3uH 

1^1 

3z 

3r  r 

'  36  J 

3z 

3v ' 

+  v|^. 

■  ♦  21 

1  w'2  ] 

3z 

3r 

r  30 

r 

3r 

aw' 

.  w' 

aw'  .  1 

.  n 

3z 

+  V- - 

3r 

+  — 
r 

ae  r 

v'w' 

-3p'  1 

ae  7 


(3) 


If  u' ,  v' ,  w' ,  p' ,  and  p'  are  each  written  as  the  sum  of  a  mean  quantity 
and  a  small  fluctuating  acoustic  quantity,  ie.,  u'  =  u  +  u,  v‘  =  v  +  v. 
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w1  *  w  ♦  w,  p'  *  p  +  p,  and  p1  ■  p  +  p,  where  the  mean  quantities  are 
Independent  of  time  and  satisfy  the  steady  mean  flow  equations  and  where 
products  of  the  fluctuating  quantities  are  considered  negligible,  then 
Equations  3  become: 

1wp  u  |jj:  ►  ~  rop 


x  vor |y.  +  + 1  aw_  +  y~|+  ap. 

YP[Jz  3r  r  36  r  J  3z  u  3r  v 

+  4  -  + vp  (|a  ♦  |3S.  +  v-  +  *w)  * 0 

36  r  H  v3z  3r  r  r  1 

(a) 

^  3u  .  w  \  ,  -  3u  .  —  3u  .  3u  u 
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In  Equations  4,  time  dependence  Is  assumed  to  be  of  the  form  e*ut  while 
9  dependence  is  assumed  to  be  of  the  form  e1me.  In  addition,  the  equation 

p 

p  «  C  p  has  been  used  to  eliminate  p  from  the  equations.  Therefore,  the 
equations  contain  only  derivatives  with  respect  to  z,  the  axial  coordinate, 
and  r,  the  radial  coordinate,  and  the  duct  can  be  taken  to  be  two- 
dimensional.  Equations  4  will  be  used  in  Section  III  in  the  finite  element 
formulation  of  the  duct  acoustics  problems  to  be  considered. 
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2.  BOUNDARY  CONDITIONS 

Because  of  the  form  of  the  differential  Equations  4,  the  appropriate 
boundary  conditions  are  the  specification  of  either  the  pressures,  the 
velocities  or  setting  the  ratio  of  the  pressure  and  velocity  equal  to  a 
given  impedance  around  the  boundary  of  the  duct.  Because  of  the  assumption 
of  axial  symmetry  along  the  part  of  the  boundary  corresponding  to  r  =  0, 
the  radial  component  of  the  velocity  is  set  to  zero.  Along  the  entrance 
of  the  duct  the  pressure  distribution  is  specified  while  at  the  exit  area 
of  the  duct  the  ratio  of  the  exit  velocity  and  pressure  are  set  equal  to 
the  exit  impedance.  These  conditions  can  be  specified  independently  of 
the  flow  within  the  duct.  However,  at  the  outer  wall  where  sound  absorbing 
liner  is  located,  the  general  form  of  the  boundary  condition  depends  on 
the  flow  if  slip  is  permitted.  That  is,  at  the  outer  wall  the  boundary 

condition  is  u  n  +  v  n  , 

o  nx  +  »  nr  •  p/z  -  i -°— C  ^  (p/t)  (5) 


where  n  =  (nx>  nr)  is  the  unit  outer  normal  to  the  boundary  at  the  point 
(x,r).  Note  that  in  the  case  of  no-slip  on  the  boundary,  uQ  nz  +  vQ  np  =  0, 
so  that  Equation  5  reduces  to  u  nx  +  v  np  =  p/z. 
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SECTION  III 

THE  FINITE  ELEMENT  METHOD 
1.  LEAST  SQUARES  METHOD 

Because  of  the  attractiveness  of  solving  large  systems  of  equations  by 
using  an  iterative  method,  the  first  attempt  at  solving  acoustic  problems 
with  a  finite  element  method  was  with  a  least  squares  method.  For  iterative 
methods  to  converge,  it  is  necessary  that  the  system  to  be  solved  be 
positive  definite  (Reference  9).  However,  if  a  straight-forward  finite 
element  approach  is  taken  with  the  acoustic  equations,  the  resulting 
system  will  not  be  positive  definite.  On  the  other  hand,  the  system  of 
equations  resulting  from  a  least  squares  approach  will  be  positive 
definite  and,  therefore,  amenable  to  an  iterative  method.  The  least 
squares  formulation  of  the  problem  consisted  of  squaring  the  left-hand 
side  of  Equations  4a-d,  summing  the  squares,  integrating  the  sum  of 
squares  over  the  area  of  the  duct,  and  then  minimizing  this  expression  as 
a  functional  of  the  variables  p,  u,  v,  w.  In  this  formulation  the  boundary 
conditions  were  imposed  by  writing  them  as  an  equation  with  vanishing 
right-hand  side,  squaring  this  equation,  integrating  it  over  the  boundary 
and  adding  this  line  integral  to  the  area  integral  already  obtained. 

The  functions  p,  u,  v,  w  were  written  as  linear  combinations  of  piece- 
wise  polynomial  functions  defined  over  rectangular  and  triangular  sub¬ 
divisions  of  the  duct,  and  these  linear  combinations  were  substituted 
into  the  aforementioned  functional.  In  this  way,  the  functional  was 
then  a  function  of  4  N  unknown  variables  instead  of  a  functional.  There 
were  4  N  variables  because  each  of  four  functions  p,  u,  v,  w  was  written 
as  a  combination  of  N  piecewise  polynomials  so  there  were  N  unknown 
coefficients  for  each  function  or  a  total  of  4  N  unknowns  to  determine. 

The  usual  procedure  of  differentiating  this  function  with  respect  to  each 
of  the  variables  and  setting  each  equal  to  zero  results  in  4  N  equations 
to  be  solved. 

Although  the  solutions  obtained  from  the  least  squares  approach  were 
qualitatively  correct  and  quantitatively  adequate,  when  compared  to  finite 
difference  solutions  requiring  approximately  the  same  computational  effort 
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they  were  not  nearly  as  accurate.  The  following  example  illustrates  the 
source  of  error  in  the  least  squares  approach  <r!nch  caused  the  poor  com¬ 
parison  with  finite  differences  when  both  were  compared  to  an  exact 
solution.  Suppose  the  ordinary  differential  equation  and  boundary  con¬ 
ditions  to  be  solved  are 

y"  +  y  =  0,  y ( 0 )  =  0,  y(l)  =  1  (6) 

on  the  interval  between  zero  and  one.  A  typical  least  squares  formulation 
would  be  to  minimize  the  functional 

J(y)  'fly"  +  y)2  dx  (7) 

o 

over  the  solution  space  of  all  functions  having  an  integrable  second 
derivative  which  vanish  at  x  =  0  and  equal  one  at  x  =  1 .  However,  if  the 
first  variation  of  this  functional  is  taken,  the  Euler  equation  obtained  is 

yIV  +  2  y1 '  +  y  =  0  (8) 

with  the  additional  natural  boundary  conditions 

(y"  +  y)  =  0  at  x=?  and  'y'  1  +  y)  =  0  at  x=l  (9) 

Now  the  general  solution  of  the  differential  Equation  6  is  y  =  a  sin  (x) 

+  b  cos  (x)  and  the  boundary  conditions  require  that  b  =  0  and  a  =  l/sin(l) 
Therefore,  the  solution  to  Equation  6  is  y  =  sin(x)/sin(l ) .  The  general 
solution  of  differential  Equation  8  is  y  =  a  sin(x)  +  b  cos(x)  +  c  x  sin(x) 
+  d  x  cos(x).  If  the  boundary  conditions,Equation  9  are  satisfied  exactly, 
then  c  =  0  and  d  =  0,  and  the  requirement  that  the  solution  vanishes  at 
x  =  0  and  equals  one  at  x  =  1  yields  y  =  sin(x)/sin  (1).  Therefore,  if 
the  natural  boundary  conditions  are  satisfied  exactly,  the  functional 
Equation  8  yields  the  correct  solution.  But  if  Equation  7  is  minimized 
approximately  using  the  finite  element  method,  then  Equation  9  will  not  be 
satisfied  exactly,  and  consequently  the  coefficients  c  and  d  will  not 
vanish.  Hence,  the  finite  element  solution  will  contain  contributions 
from  the  two  spurious  functions  y3  =  x  sin(x)  and  y4  =  x  cos(x).  This  is 
the  additional  error  that  one  encounters  when  using  a  least  squares 
finite  element  method  which  is  basic  to  the  problem  and  not  caused  by 
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nonconforming  elements.  As  a  result  of  this  additional  error,  a  finer 
discretization  is  required  than  would  be  necessary  if  it  were  not  present. 

If  the  iterative  solution  of  the  least  squares  problem  were  sufficiently 
quicker  in  terms  of  computer  time  to  offset  the  additional  computational 
work  required  as  a  result  of  the  finer  mesh,  then  the  least  squares  method 
would  still  be  preferable  to  a  Galerkin  method  which  requires  a  direct 
system  solver.  To  compare  the  two  possibilities,  the  equations  were  solved 
using  both  the  iterative  approach  and  direct  Gaussian  elimination  and  the 
computational  times  were  compared.  For  a  system  with  314  unknowns  and  a 
bandwidth  of  42  (42  non-zero  entries  in  any  given  row  and  all  of  these 
entries  clustered  about  the  main  diagonal),  a  moderately  accurate  iterative 
solution  was  obtained  with  41  iterations  taking  14  seconds  of  CP  time 
and  373  seconds  of  input-output  time  on  the  compute''.  For  the  same  problem, 
the  direct  Gaussian  elimination  took  only  3  seconds  of  CP  time  and  10 
seconds  of  input/output  time.  The  conclusion  was,  therefore,  that  for 
the  size  matrices  considered,  a  direct  solver  was  preferred  to  an  iterative 
method.  Because  the  major  attraction  of  the  least  squares  method  was  the 
possibility  of  solving  the  large  system  of  equations  iteratively  and  since 
the  iterative  methods  did  not  turn  out  to  be  less  expensive  computationally, 
it  was  decided  to  use  a  Galerkin  or  weak  solution  method.  The  integral  or 
weak  form  of  the  equations  will  now  be  described. 

2.  THE  INTEGRAL  FORM  OF  THE  EQUATIONS 


To  avoid  rewriting  Equations  4  repeatedly,  it  will  be  represented  by: 

Dl(p,u,v,w)  =0  (a) 

D2(p,u,v,w)  =  0  (b) 

(4) 

D3(p,u,v,w)  =  0  (c) 

D4(p,u, v,w)  =  0  (d) 


Then  the  weak  or  integral  form  of  Equations  4  is  obtained  by  multiplying 
each  of  the  equations  by  the  product  of  the  variable  r  and  a  test  function 
(weight  function)  and  integrating  over  the  area  of  the  duct.  In  this  way 
the  requirement  that  the  solutions  of  Equations  4  have  continuous  first 
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derivatives  within  the  duct  is  relaxed  to  the  requirement  that  they  have 
generalized  first  derivatives  which  are  integrable.  The  integral  form  of 
Equations  4  is  written  as 

jj  r  f  Dl(p,u,v,w)  dz  dr  =  0  U) 

//  r  f  D2(p,u,v,w)  dz  dr  =  0  (b)  ^ 

//  r  f  D3(p,u,v,w)  dz  dr  =  0  (c) 

//  r  f  D4(p,u,v,w)  dz  dr  =  0  (d) 

where  f  is  a  test  function.  The  boundary  conditions  relating  the  ratio 
of  the  pressure  and  velocity  to  the  impedance  are  specified  by  integrating 
some  of  the  terms  in  Equation  10a  by  parts. 

The  particular  integration  by  parts  is 

//  Yf  P0  (r  •—  +  r  |^-  +  v  +  m  w)  dz  dr 

=  /yr  f  pQ  (u  nz  +  v  nr)  ds  +  //jyf  pQ  m  w  (11) 

-  yr  <|fo  f  +  p0  |i  )u-yr  (fpo  f  +  P0  If)']  dz  dr 

Therefore,  the  term  in  Equation  10a  which  corresponds  to  the  left-hand 
side  of  Equation  11  is  replaced  by  the  right-hand  side  of  Equation  11. 

This  new  version  of  Equation  10a  will  be  denoted  as  Equation  10a'. 

A  finite  element  solution  of  Equations  10a',  b,  c,  and  d  consists 
of  first  dividing  the  interior  of  the  duct  into  a  collection  of  rectangles 
and  triangles  and  then  distributing  an  appropriate  number  of  grid  points 
or  nodes  on  the  sides  and  in  the  interior  of  these  triangles  and 
rectangles.  Figure  1  indicates  a  typical  finite  element  discretization 
of  two  difference  duct  shapes.  If  piecewise  linear  basis  functions  (shape 
functions)  are  to  be  used  on  triangles  then  a  node  at  each  of  the  vertices 
of  each  of  the  triangles  is  appropriate.  If  piecewise  linear  basis 
functions  are  to  be  used  on  rectangles  then  a  node  at  each  of  the  vertices 
of  each  rectangle  is  appropriate.  A  node  at  each  vertex,  a  node  at  the 
midpoint  of  each  side  and  a  node  at  the  center  of  each  rectangle  for  a 
total  of  nine  nodes  is  appropriate  if  piecewise  biquadratic  basis  functions 
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are  to  be  used.  In  Figure  2,  the  location  of  nodes  on  triangles  and  rec¬ 
tangles  is  indicated.  Whatever  the  choice  of  basis  function,  the  next  step 
Is  to  write  p,u,v,w  each  as  a  linear  combination  of  the  basis  functions. 
That  is. 
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where  each  g.  is  a  basis  function,  the  p.,  u.,  v.  and  w.  are  unknown 
J  J  J  «  J 

coefficients  to  be  determined,  and  N  is  the  number  of  basis  functions. 


Figure  la.  Discretization  of  a  Conical  Duct 

11 


»■  •  ■  * '  — r  1 1 1  if ii «■<  I  t friV *i >  ■ a ^  - 


AFWAL-TR-81-3087 


The  are  determined  by  requiring  the  jth  function  to  be  1  at  the  jth 

V 

node  and  to  vanish  at  all  other  nodes.  Typically  this  is  done  on  an 
element  basis;  that  is,  the  nine  biquadratic  basis  functions  which  are 
nonzero  within  a  rectangle  are  determined  by  specifying  each  to  be  one  at 
one  node  on  the  rectangle  and  vanish  at  the  other  eight  nodes.  In  this 
way  a  system  of  equations  Is  obtained  which  determines  the  coefficients 
of  the  nine  biquadratic  basis  functions.  Once  the  g.  have  been  deter- 

J 

mined,  the  expression  for  p,u,v,w  dtfined  by  Equations  12  are  substituted 
into  Equations  10a',  b,  c,  and  d.  The  test  function  f  which  appears  in 
Equations  10  is  set  equal  to  each  of  the  functions  g.  for  j  *  1  through 

J 

j  =  N.  In  this  way  4  N  algebraic  equations  are  obtained  which  determine 
the  4  N  unknown  coefficients  p^ ,  ...,  p^,  u^ ,  u^,  v-j ,  ...,  vN, 

w^ ,  ....  wN.  Then  if  the  value  of  p,u,v  or  w  is  needed  anywhere  within 
the  duct.  Equation  12  is  used  with  the  computed  values  of  the  coefficients 
in  Equation  12. 

Because  of  the  large  system  of  equations  which  arises  in  the  deter¬ 
mination  of  the  p.'s,  u.'s,  v.'s,  and  w.'s,  an  out-of-core  solver  is 

J  J  J  J 

almost  mandatory.  Even  for  a  relatively  coarse  discretization  of  5  x  5 
elements  with  biquadratic  basis  functions,  there  are  121  nodes  and  con¬ 
sequently  484  equations  to  be  solved.  Even  with  a  system  solver  which 
accounts  for  a  bandwidth  of  only  200,  to  keep  the  entire  banded  system, 
right-hand  side  and  auxiliary  vectors  in  core  requires  96,800  memory 
locations.  This  is  almost  300,000  octal  locations,  the  maximum  core 
available  with  the  CDC  CYBER  machines  at  the  ASD  Computer  Center  at 
Wright-Patterson.  Therefore,  an  out-of-core  banded  solver  is  used  which 
requires  a  moderate  amount  of  in-core  storage  but  much  more  input/output 
time. 

Because  of  the  ease  of  switching  from  linear  to  quadratic  basis 
functions  (Figure  3  shows  the  ease  In  changing  fm:r.  one  mesh  to  the  other), 
it  was  decided  to  compare  the  accuracy  of  the  two  types  of  basis  functions 
in  deciding  which  to  use.  The  following  problem  was  solved  on  the  same 
grid  on  a  rectangular  duct.  The  differential  equation  was  the  Helmholtz 
equation  written  in  system  form. 
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Figure  4.  Least  Squares  Error  as  a  Function  of  Grid  Size  for  Linear 
and  Quadratic  Basis  Functions 


.....  9pl  3pl  3P2  3p3  2 

chat  is,  p2  =  P3  ■  +  a^Pl  -  0,  with  p,  =  1  for 

x  =  0,  p3  =  0  tor  v  *  0  and  for  y  =  1 ,  and  p2  +  i  a  p1  =  0  for  x  *  1 . 
The  solutions  obtained  uy  using  the  bilinear  basis  functions  and  the 
biquadratic  basis  function:,  were  compared  to  the  exact  solution  for 
several  different  grids,  and  the  error  is  plotted  as  a  function  of  grid 
size  in  Figure  4.  The  error  for  bilinear  basis  functions  varies  as  the 
second  power  of  the  mesh  size,  while  the  error  for  the  biquadratic  basis 
functions  vanes  as  the  fourth  puwer  of  the  mesh  size.  Therefore, 
reducing  the  mesh  size  hy  a  factor  ot  two  reduces  the  error  by  a  factor 
of  four  if  bilinear  basis  functions  are  used  while  it  reduces  the  error 
by  a  factor  of  sixteen  when  biquadratic  basis  functions  are  used.  As  a 
result  of  this  experience,  it  was  decided  to  use  the  biquadratics  and 
obtain  the  improved  accuracy  with  less  computational  time. 
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One  final  decision  that  was  made  was  how  to  handle  the  nonuni  form 
geometries.  The  most  straightforward  approach  was  to  use  rectangles  and 
triangles  with  the  triangles  used  to  approximate  the  boundary  as  indicated 
in  Figure  1.  However,  as  a  test,  the  case  of  rto-flow  in  a  hard-walled 
uniform  duct  was  considered.  With  a  plane  wave  as  input  and  an  appropriate 
exit  Impedance,  the  exact  solution  would  be  p(z,r)  *  cos(az)-1  sin(az) 
where  a  is  the  specified  frequency.  With  rectangular  elements,  the  finite 
element  solution  duplicated  the  exact  solution,  particularly  in  the  respect 
that  the  solution  was  independent  of  the  variable  r;  that  is,  there  was  no 
variation  in  the  r  direction.  When  triangular  elements  were  specified  by 
dividing  each  rectangle  into  two  triangles  as  indicated  in  Figure  5,  the 
error  Increased,  and  the  finite  element  solution  varied  incorrectly  as  a 
function  of  r.  This  is  consistent  with  the  observations  of  Reference  4. 
With  this  experience,  it  was  decided  to  look  closely  at  the  isoparametric 
transformations  and  the  use  of  quadrilaterals  with  curved  boundaries. 

In  Figure  6,  the  two  ducts  which  were  depicted  in  Figure  1  are  subdivided 
into  curve-sided  quadrilaterals.  The  two  main  difficulties  with  using 
these  elements  is  the  definition  of  the  biquadratics  on  such  a  nonuniform 
geometry  and  the  integration  over  such  a  geometry.  The  advantages  are 
the  Improved  approximation  of  curved  boundaries  and  the  ease  of  generating 
the  elements  or  subdivisions.  Th.'  two  disadvantages  are  easily  overcome 
by  mapping  the  nonuniform  element  to  a  square  as  indicated  in  Figure  7. 

On  the  square  the  mapped  biquadratics  become  the  usual  biquadratics,  and 
the  integration  is  straightforward.  The  apparent  difficulty  of  deter¬ 
mining  the  mapping  function  is  also  easily  resolved.  The  step  which 
makes  the  mapping  easy  is  the  approximation  of  the  original  nonuniform 
geometry  by  piecewise  quadratics  around  the  boundary  and  requiring  the 
sides  of  the  quadrilateral  elements  with  curved  sides  to  be  quadratic 
curves.  Then  the  map  from  the  (r, ,n)  plane  to  the  (z,r)  plane  in  Figure  7 
is  the  biquadratic  which  maps  each  of  the  nine  nodes  in  the  U,n)  plane 
to  the  appropriate  node  of  the  nonuniform  element  in  the  (z,r)  plane. 

But  any  biquadratic  on  the  U,n)  square  can  be  written  as  a  linear 
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combination  of  the  nine  elementary  biquadratic  basis  functions  defined 
therr.  That  Is,  If  the  map  from  the  (c,n)  plane  to  the  (z.r)  plane  Is 
z  «  z  (c,n)  and  r  «  r  (s,n)  then  i  and  r  can  be  written  as: 

9 

z*z(S,nH  a,  g*U.n) 

,1-1  3  3 

9 

r-rU.n)*!  b,  g.U.n)  (13) 

j*l  3  3 


To  determine  the  nine  coefficients  a.,  j  «  1  to  j  *  9,  and  the  nine 

J 

coefficients  bj,  j  *  1  to  j  *  9,  one  needs  to  recall  that  zk  *  *Uk*nk) 


and  that 


3  £,  aj  9j  th,t  rk  *  r(W3£,  bj9j«ek*nk) 

gjUk,nk)  a  1  if  jak»  and  9jUk.nk)  a  0  If  jfk.  Therefore  a^  ■  z.  and 

bj  ■  Tj,  and  the  map  Is  determined, 

A  typical  Integration  In  the  original  (z.r)  plane  Is  transformed 
Into  an  Integration  In  the  U.n)  plane  as  follows: 


‘  t 

. 


f (z.r)  p(z,r)  dz  dr 
11 


3fl< 

00 


f(z(E;,n).rU,n))  p(zU.n)  »r(E^ ) 


(14) 


where  J  is  the  Jacobian  of  the  transformation  and  is  given  by 
i  -  12.  it  i2  it. 

J  "  35  3n  "  3n  35  •  The  partial  derivatives  of  c  and  n  on  the  right- 
hand  side  of  Equation  14  are  given  by: 


3£  _  3r.,  3n  *  -3r «  ,  3E  _  -3Z/,  ,_j 

9,  ‘  9?  3  If  -  S/J'  a"d 


9n.  M.  ,, 
3r  H  /J* 


The  actual  Integration  on  the  right-hand  side  of  Equation  14  Is  performed 
with  three-point  Gaussian  quadrature  In  each  of  the  coordinate  directions 
for  a  total  of  nine  quadrature  points  within  the  square.  The  results  for 
uniform  and  nonuniform  ducts  with  and  without  flow  will  now  be  discussed. 
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SECTION  IV 
RESULTS 

The  first  test  case  for  the  finite  element  program  was  to  compute 
propagation  of  a  plane  wave  in  a  hard-walled  uniform  duct  not  containing 
flow.  The  outstanding  agreement  with  the  exact  solution  is  indicated  in 
Figure  8.  For  a  second  test  case  the  variation  of  an  individual  mode  in 
a  soft  wall  duct  not  containing  flow  was  computed.  As  seen  in  Figure  9, 
the  agreement  between  the  finite  element  solution  and  the  exact  solution 
is  good  but  not  as  good  as  for  the  plane  wave.  This  is  due  to  the 
variation  in  the  radial  direction.  For  the  cases  of  uniform  flow  in  a 
uniform  duct,  a  Mach  number  of  .2  with  both  a  plane  wave  in  a  hard  wall 
duct  and  a  single  mode  in  a  soft-walled  duct  was  considered.  In  Figures 
10  and  11  the  comparisons  between  the  finite  element  and  exact  solutions 
are  indicated.  Again,  the  agreement  is  good,  with  agreement  for  the  plane 
wave  in  the  hard  wall  duct  outstanding.  For  a  nonuniform  duct  the  two 
geometries  of  Figure  6  were  considered  to  be  typical.  The  cone  was  most 
interesting  because  the  exact  solution  for  a  plane  wave  propagating  in  a 
hard-walled  cone  is  known  and  because  the  finite  difference  conformal  map 
solution  is  available.  In  Figures  12a  and  12b  the  results  for  a  plane  wave 
propagating  in  a  hard-walled  cone  for  two  different  discretizations  are 
depicted.  For  the  cone,  the  agreement  with  exact  solution  is  not  as  good 
for  the  finite  element  method  as  for  the  finite  difference  method.  This 
is  because  the  finite  difference  method  uses  an  exact  map  to  a  uniform 
duct  while  the  finite  element  solution  method  approximates  the  circular 
inlet  and  exit  as  piecewise  quadratic.  This  is  shown  by  the  decrease  in 
accuracy  in  Figure  12b  where  the  number  of  elements  in  the  radial 
direction  was  reduced  in  comparison  with  the  number  in  Figure  12a. 
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